Why Use R?
“R is a free software environment for statistical computing and graphics”
Compared to Commercial Flow Cytometry Software, R has the following advantages:- Reproducible (Data + Code = Figures)
- Flexible (This presentation was created in R!)
- Nice Graphics (ggplot2)
- Great for analysis pipelines and frequently used workflows
- Newest analysis techniques
- Free!
Representing Data in R
- Matrix: An \(n*m\) array of items, all of the single class
- Data Frame: An \(n*m\) array of items, but each column can be a different class
Example: The iris dataset: measurements of 50 flowers of 3 species of iris
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
Defining Tidy Data
To work with data in R, it’s best to have ‘tidy data,’ which meets the following criteria:
- Each variable [‘feature’] must have its own column.
- Each observation [‘cell’] must have its own row.
- Each value must have its own cell [‘entry’].
…but cytometry data is not usually tidy.
For more information, see: Wickham, Hadley. “Tidy data.” Journal of Statistical Software 59.10 (2014): 1-23.
Making Cytometry Figures in R (1)
One thing we may want to do is reproduce the same tSNE figure we made on cytobank:
ggplot(mydataframe, aes(x = tSNE1, y = tSNE2)) +
geom_point(shape = ".") +
coord_fixed() +
facet_wrap(~Group) 
Making Cytometry Figures in R (2)
We can also customize our plots in ways that are not easy to do in cytobank:
ggplot(mydataframe, aes(x = tSNE1, y = tSNE2)) +
geom_point(shape = 16, alpha = 0.2, size = 0.2) +
coord_fixed() +
facet_wrap(~Group) +
scale_color_viridis_c(option = "A") +
theme_minimal() +
theme(axis.text = element_blank())
Making Cytometry Figures in R (4)
Then we’ll make our plot:
mydataframe.long %>%
ggplot(aes(x = tSNE1, y = tSNE2, col = intensity)) +
geom_point(shape = ".") +
scale_x_continuous(expand = c(0,0)) +
scale_y_continuous(expand = c(0,0)) +
scale_color_viridis_c(option = "A") +
coord_fixed() +
facet_wrap(~marker, nrow = 3) +
theme_minimal() +
theme(axis.text = element_blank())
Applying alternative dimensionality reduction techniques (3)
Next, we’ll rejoin the two new UMAP columns to our original dataframe, and make our plot:
Applying alternative dimensionality reduction techniques (4)
mydataframe %>%
bind_cols(as.data.frame(myumap)) %>%
as_tibble() %>%
gather(marker, intensity, contains("(V)")) %>% # <- this is the key step
separate(marker, c("channel", "marker", "drop"), sep= "_|\ ") %>%
ggplot(aes(x = UMAP1, y = UMAP2, col = intensity)) +
geom_point(shape = ".") +
scale_x_continuous(expand = c(0,0)) +
scale_y_continuous(expand = c(0,0)) +
scale_color_viridis_c(option = "A") +
coord_fixed() +
facet_wrap(~marker, nrow = 3) +
theme_minimal() +
theme(axis.text = element_blank())
Applying alternative dimensionality reduction techniques (5)
We can also plot our data as a map:
library(scico)
axis.max <- apply(myumap, 2, max) + 1
axis.min <- apply(myumap, 2, min) - 1
mydataframe %>%
bind_cols(as.data.frame(myumap)) %>%
ggplot(aes(x=UMAP1, y = UMAP2)) +
stat_density_2d(h = c(1, 1),
n = 1024,
geom = "raster",
contour = F,
aes(fill = stat(density))) +
scale_fill_scico(palette = "oleron", name = "density", trans = "sqrt") +
scale_x_continuous(expand = c(0,0), limits = c(axis.min[1], axis.max[1])) +
scale_y_continuous(expand = c(0,0), limits = c(axis.min[2], axis.max[2])) +
theme_classic() +
coord_fixed()
Clustering in R
We can also apply a clustering algorithm to our dimensionality reduced data.
- Density-based spatial clustering of applications with noise (DBSCAN)
- No need to specify number of clusters
- Few parameters to tune (eps and minPts)
- Fast + scalable
#install.packages("dbscan")
library(dbscan)
mydbscan <- dbscan(myumap, eps = 0.3, minPts = 150)
mydbscan## DBSCAN clustering for 100000 objects.
## Parameters: eps = 0.3, minPts = 150
## The clustering contains 22 cluster(s) and 453 noise points.
##
## 0 1 2 3 4 5 6 7 8 9 10 11
## 453 2073 8279 37863 4422 6956 2559 15082 1102 1219 2836 2133
## 12 13 14 15 16 17 18 19 20 21 22
## 2693 6472 376 1731 1722 677 522 208 353 176 93
##
## Available fields: cluster, eps, minPtsmydataframe %>%
bind_cols(as.data.frame(myumap)) %>%
mutate(cluster = as.factor(mydbscan$cluster)) %>%
# dplyr::filter(cluster == 1:3) %>%
ggplot(aes(x=UMAP1, y = UMAP2, col = cluster)) +
geom_point(shape = ".") +
scale_color_discrete(guide = guide_legend(override.aes = list(shape = 15))) +
coord_fixed() +
theme_classic()
Understanding our clusters
We have clusters, but how can we understand what makes them distinct?
- Marker Enrichment Modeling (MEM)
- Heatmaps (see example below)
library(tibble)
library(scico)
library(seriation)
myheatmap <- mydataframe %>%
bind_cols(as.data.frame(myumap)) %>%
mutate(cluster = as.factor(mydbscan$cluster)) %>%
gather(marker, intensity, contains("(V)")) %>% # <- this is the key step
group_by(cluster, marker) %>%
summarise(MFI = median(intensity)) %>%
select(marker, MFI, cluster) %>%
spread(marker, MFI)
myheatmap.mat <- myheatmap %>%
ungroup() %>%
column_to_rownames("cluster") %>%
as.matrix()
matrix.dist.row <- dist((myheatmap.mat))
matrix.dist.col <- dist(t(myheatmap.mat))
row.order <- seriation::get_order(seriate(matrix.dist.row, method = "HC"))
col.order <- seriation::get_order(seriate(matrix.dist.col, method = "HC"))
myheatmap %>%
ungroup() %>%
gather(marker, MFI, -cluster) %>%
mutate(marker = factor(marker, levels = colnames(myheatmap.mat)[(col.order)])) %>%
mutate(cluster = factor(cluster, levels = rownames(myheatmap.mat)[(row.order)])) %>%
ggplot(aes(x=cluster, y = marker,fill = MFI)) +
geom_tile(color = "grey90", size = 0.5) +
scale_fill_viridis_c(option = "E") +
scale_x_discrete(expand = c(0,0)) +
scale_y_discrete(expand = c(0,0)) +
coord_fixed()
What else is out there
- Chen, H, et al., PLOS Computational Biology 2016
- Integrated pipeline for analysing cytometry data in R
- Set of interactive tools that integrate with flowWorkspace
- Interactive gating, compensations, panels, etc…
